#include "LinearSystem.h"

/*
 *--------------------------------------------------------------------------------------
 *       Class:  LinearSystem
 *      Method:  LinearSystem
 * Description:  Contrói um sistema linear vazio
 *--------------------------------------------------------------------------------------
 */
LinearSystem::LinearSystem()
{
    order = testRow = error = maxIt = 0;
    A = b = x0 = x1 = NULL;
}

/*
 *--------------------------------------------------------------------------------------
 *       Class:  LinearSystem
 *      Method:  ~LinearSystem
 * Description:  Destrói o sistema linear
 *--------------------------------------------------------------------------------------
 */
LinearSystem::~LinearSystem()
{
    if(A != NULL)
        delete A;
    if(b != NULL)
        delete b;
    if(x0 != NULL)
        delete x0;
    if(x1 != NULL)
        delete x1;
}

/*
 *--------------------------------------------------------------------------------------
 *       Class:  LinearSystem
 *      Method:  solve
 * Description:  Resolve o sistema linear pelo método iterativo de Jacobi
 *--------------------------------------------------------------------------------------
 */
void LinearSystem::solve()
{
    mdouble& A = REF(A);
    mdouble& b = REF(b);
    mdouble& x0 = REF(x0);
    mdouble& x1 = REF(x1);

    bool stop = false;

    /* Loop que para quando há a convergência ou quando atinge a iteração máxima */
    for(iterations = 0; !stop && iterations < maxIt; iterations++)
    {
        /* Faz a conta x1(i) = 1/A(i,i) * [b(i) - sum{A(i,j) * x0(j)}] */
        for(uint i = 0; i < order; i++)
        {
            x1(i) = 0;

            for(uint j = 0; j < order; j++)
                x1(i) += A(i, j) * x0(j);

            x1(i) = x0(i) + (b(i) - x1(i)) / A(i, i);
        }

        /* Verifica se convergiu */
        stop = hasConverged();

        /* Shifta os valores das variáveis */
        for(uint i = 0; i < order; i++)
            x0(i) = x1(i);
    }
}

/*
 *--------------------------------------------------------------------------------------
 *       Class:  LinearSystem
 *      Method:  parallelSolve
 * Description:  Resolve o sistema linear pelo método iterativo de Jacobi de forma paralela
 *--------------------------------------------------------------------------------------
 */
void LinearSystem::parallelSolve()
{
    mdouble& A = REF(A);
    mdouble& b = REF(b);
    mdouble& x0 = REF(x0);
    mdouble& x1 = REF(x1);

    /* Número de processos que estão fazendo os cálculos */
    uint nProcs = Manager::getNProcs() - 1;

    char stop = 0;

    /* Loop que para quando há a convergência ou quando atinge a iteração máxima */
    for(iterations = 0; !stop && iterations < maxIt; iterations++)
    {
        /* Recupera o número de linhas e a linha de início desse processo */
        uint nRows = Manager::getNRows();
        uint initRow = Manager::getInitRow();

        for(uint step = 0; step < nProcs; step++)
        {
            /* Recupera o número de linhas e a linha de início nesse passo */
            uint currNRows = Manager::getNRows(step);
            uint currInitRow = Manager::getInitRow(step);

            /* Para cada x1 responsável */
            #pragma omp parallel for
            for(uint i = initRow; i < initRow + nRows; i++)
                /* Para cada x0 sob posse atualmente */
                for(uint j = currInitRow; j < currInitRow + currNRows; j++)
                    x1(i) += A(i, j) * x0(j);

            /* Envia os valores de x0 para o processo anterior */
            x0.send(Manager::prevRank(), currInitRow, currNRows);

            /* Recebe os valores de x0 do processo posterior e espera o recebimento */
            x0.recv(Manager::nextRank(), Manager::getInitRow(step+1), Manager::getNRows(step+1));
        }

        /* Finaliza a fórmula de x1 e shifta os valores de x1 para x0 */
        #pragma omp parallel for
        for(uint i = initRow; i < initRow + nRows; i++)
        {
            x0(i) = x0(i) + (b(i) - x1(i)) / A(i, i);
            x1(i) = 0;
        }

        /* Envia os valores calculados para o nó mestre */
        x0.send(0, initRow, nRows);

        /* Verifica se convergiu (obtém valor do no mestre) */
        recvStopSignal(&stop, 0);
    }
}

/*
 *--------------------------------------------------------------------------------------
 *       Class:  LinearSystem
 *      Method:  parallelCheck
 * Description:  Verifica por convergência paralelamente pelo nó mestre
 *--------------------------------------------------------------------------------------
 */
void LinearSystem::parallelCheck()
{
    mdouble& x0 = REF(x0);
    mdouble& x1 = REF(x1);

    /* Número de processos que estão fazendo os cálculos */
    uint nProcs = Manager::getNProcs();

    char stop = 0, warned = 0;

    for(iterations = 0; !warned && iterations < maxIt; iterations++)
    {
        /* Recupera todos os valores calculados */
        for(uint i = 1; i < nProcs; i++)
        {
            sendStopSignal(&stop, i); /* Envia sinal de stop aos escravos */
            x1.recv(i, Manager::getInitRow(i), Manager::getNRows(i));
        }
        warned = stop;

        if(warned)
            continue;

        /* Verifica se terminou o algoritmo */
        if(iterations >= 1)
            stop = hasConverged();

        /* Shifta os valores de x1 para x0 */
        #pragma omp parallel for
        for(uint i = 0; i < order; i++)
            x0(i) = x1(i);
    }

    iterations--;
}

/*
 *--------------------------------------------------------------------------------------
 *       Class:  LinearSystem
 *      Method:  hasConverged
 * Description:  Verifica se o método de Jacobi já convergiu
 *--------------------------------------------------------------------------------------
 */
//bool LinearSystem::hasConverged()
//{
//    mdouble& x0 = REF(x0);
//    mdouble& x1 = REF(x1);
//
//    double max = fabs(x1(0) - x0(0));
//
//    /* Pega a maior diferença entre os dois vetores */
//    for(uint i = 1; i < order; i++)
//    {
//        double d = fabs(x1(i) - x0(i));
//        if(d > max)
//            max = d;
//    }
//
//    /* Verifica se ele é menor que o erro máximo permitido */
//    return (max < error);
//}

bool LinearSystem::hasConverged()
{
    mdouble& b = REF(b);

    double diff = fabs(b(testRow) - calcRow());

    /* Verifica se ele é menor que o erro máximo permitido */
    return (diff < error);
}

/*
 *--------------------------------------------------------------------------------------
 *       Class:  LinearSystem
 *      Method:  sendStopSignal
 * Description:  Envia sinal de parada (convergência) do algoritmo
 *--------------------------------------------------------------------------------------
 */
void LinearSystem::sendStopSignal(char *stop, uint to)
{
    MPI::COMM_WORLD.Isend(stop, 1, MPI::CHAR, to, 2);
}

/*
 *--------------------------------------------------------------------------------------
 *       Class:  LinearSystem
 *      Method:  recvStopSignal
 * Description:  Recebe sinal de parada (convergência) do algoritmo
 *--------------------------------------------------------------------------------------
 */
void LinearSystem::recvStopSignal(char *stop, uint from)
{
    MPI::COMM_WORLD.Recv(stop, 1, MPI::CHAR, from, 2);
}

/*
 *--------------------------------------------------------------------------------------
 *       Class:  LinearSystem
 *      Method:  calcRow
 * Description:  Testa a linha desejada
 *--------------------------------------------------------------------------------------
 */
double LinearSystem::calcRow()
{
    mdouble& A = REF(A);
    mdouble& x1 = REF(x1);

    double sum = 0;

    /* Calcula o valor da linha de teste */
    for(uint j = 0; j < order; j++)
        sum += A(testRow, j) * x1(j);

    return sum;
}

/*
 *--------------------------------------------------------------------------------------
 *       Class:  LinearSystem
 *      Method:  operator >>
 * Description:  Adiciona os valores do sistema a partir de um stream de entrada
 *--------------------------------------------------------------------------------------
 */
istream& operator >>(istream& is, LinearSystem& ls)
{
    is >> ls.order >> ls.testRow >> ls.error >> ls.maxIt;

    ls.A = new mdouble(ls.order, ls.order);
    ls.b = new mdouble(ls.order);
    ls.x0 = new mdouble(ls.order);
    ls.x1 = new mdouble(ls.order);

    /* Iniciliza o vetor inicial com 0 */
    ls.x0->fillZeros();

    is >> *(ls.A) >> *(ls.b);

    return is;
}

/*
 *--------------------------------------------------------------------------------------
 *       Class:  LinearSystem
 *      Method:  operator <<
 * Description:  Imprime os valores em um stream de saída
 *--------------------------------------------------------------------------------------
 */
ostream& operator <<(ostream& os, LinearSystem& ls)
{
    os << "---------------------------------------------------------" << endl;
    os << "Iterations: " << ls.iterations << endl;
    os << "RowTest: " << ls.testRow << " => [" << ls.calcRow() << "] =? " << (*ls.b)(ls.testRow) << endl;
    os << "---------------------------------------------------------";

    return os;
}
